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We present numerical and analytical results for the optical and acoustic plasmon collective modes 
of coupled massless-Dirac two-dimensional electron systems. Our results apply to topological insu- 
lator (TI) thin films and to two graphene sheets separated by a thin dielectric barrier layer. We 
find that because of strong bulk dielectric screening TI acoustic modes are locked to the top of the 
particle-hole continuum and therefore probably unobservable. 
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I. INTRODUCTION 

The physics of closely-spaced but unhybridized two- 
dimensional electron systems (2DESs) has been a sub- 
ject of theoretical and experimental interest since it was 
first appreciated^ ^ that electron-electron interactions al- 
low energy and momentum to be transferred between 
layers, while maintaining separate particle-number con- 
servation. Remote Coulomb coupling has commanded a 
great deal of attention during the past thirty years or so 
because it provides a potential alternative to the induc- 
tive and capacitive coupling of conventional electronics. 
Until recently, remote Coulomb coupling research focused 
on quasi-2D electron systems confined to nearby quan- 
tum wells in molecular-beam-epitaxy grown semiconduc- 
tor heterostructures. The study of Coulomb-coupled 2D 
systems has now been revitalized by advances which have 
made it possible to prepare robust and ambipolar 2DESs, 
based on graphene^ layers or on the surface states of 
topological insulators^, that are described by an ultra- 
relativistic wave equation instead of the non-relativistic 
Schrodinger equation. 

Single- and few-layer graphene systems can be pro- 
duced by mechanical exfoliation of thin graphite or 
by thermal decomposition of silicon carbide^. Isolated 
graphene layers host massless-Dirac two-dimensional 
electron systems (MD2DESs) with a four-fold (spin 
X valley) flavor degeneracy, whereas topologically- 
protected MD2DESs that have no additional spin or val- 
ley flavor labels appear automaticall}BEl ^t the top and 
bottom surfaces of a three-dimensional (3D) TI thin film. 
The protected surface states of 3D TIs are associated 
with spin-orbit interaction driven bulk band inversions. 
3D TIs in a slab geometry offer two surface states that 
can be far enough apart to make single-electron tun- 
neling negligible, but close enough for Coulomb inter- 
actions between surfaces to be important. Unhybridized 
MD2DES pairs can be realized in graphene by separat- 
ing two layers by a dielectric!^ (such as AI2O3) or by a 
few layers of a one-atom-thick insulator such as BN^. 
In both cases inter-layer hybridization is negligible and 
the nearby graphene layers are, from the point of view of 



single-particle physics, isolated. Isolated graphene lay- 
ers can be also found on the surface of bulk graphite^^'^ 
and in "folded graphene"'!^ (a natural byproduct of mi- 
cromechanical exfoliation), or prepared by chemical va- 
por deposition^^. We use the term double-layer graphene 
(DLG) to refer to a system with two graphene layers that 
are coupled only by Coulomb interactions, avoiding the 
term bilayer graphene which typically refers to two ad- 
jacent graphene layers in the crystalline Bernal-stacking 
configuration^^. 

DLG and TI thin films are both described at low ener- 
gies by a Hamiltonian with two MD2DES^ coupled only 
by Coulomb interactions. The importance of electron- 
electron interactions in MD2DESs has been becoming 
more obvious as sample quality has improved^^, motivat- 
ing investigations of charge and spin or pseudospin dy- 
namics in DLG and thin-film TIs in the regime in which 
long-range Coul omb forces give rise to robust plasmon 
collective mode4^^^. Because of their electrically tun- 
able collective behaviors, DLG and thin-film TIs may 
have a large impact on plasmonics, a very active subfield 
of optoelectronicJI^KI^ whose aim is to exploit plasmon 
properties in order to compress infrared electromagnetic 
waves to the nanometer scale of modern electronic de- 
vices. 

In this Article we use the random phase approximation 
(RPA)^^ to evaluate the optical and acoustic plasmon 
mode dispersions in DLG and in thin-film TIs. In partic- 
ular, we obtain an exact analytical formula for the RPA 
acoustic plasmon group velocity valid for arbitrary sub- 
strate and barrier dielectrics that points to a key dif- 
ference between these two MD2DES's, namely that the 
velocity in TI thin films is strongly suppressed. The RPA 
collective modes of DLG have been calculated earlier by 
Hwang and Das Sarma^^: below we will comment at 
length on the relation between our results and theirs. 
Plasmon collective modes formed from TI surface states 
have also been considered previously by Raghu et al^ 
in the regime in which coupling between top and bottom 
surfaces can be neglected. Based on our analysis, we are 
able to clarify how dielectric screening influences plasma 
frequencies in this limit. 

Plasmons can be observed by a variety of experi- 
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FIG. 1: (Color online) A side view of the double- layer system 
described by Eq. ([T]), which explicitly indicates the dielec- 
tric model used in these calculations. The two layers hosting 
massless Dirac fermions are located at z — and z — d. 



a = {(T^^cjy) is a vector of Pauli matrices. A sum over 
flavor labels is implicit in Eq. ([2| in the case of DLG. The 
relative strength of Coulomb interactions is measured by 
the dimensionless coupling constant^ (restoring h for a 
moment) = /{hv) which has a value ~ 2.2 in DLG 
and ^ 4.4 in Bi2Te3 TIs if we use the respective Dirac 
velocities vq ~ 10^ m/s and v^i 5 x 10^ m/s. 

Several important many-body properties of the Hamil- 
tonian 1-L are completely determined by the 2x2 symmet- 
ric matrix %(<7, cj) whose elements are the density-density 
linear-response functions 



mental tools including inelastic light scattering^^, which 
has been widely used to probe plasmons in semicon- 
ductor heterostructure^^, but also by surface-physics 
techniques like high-resolution electron-energy-loss spec- 
troscop}'^^, and, more indirectly, angle-resolved photoe- 
mission spectroscop}'^^. Double-layer field-effect tran- 
sistors with a grating gat^^ can also be used to de- 
tect plasmons. Coupling between far-infrared light and 
Dirac plasmons in single-layer graphene has recently been 
achieved by employing an array of graphene nanorib- 
bonP^ and by performing near-field scanning optical mi- 
croscopy through the tip of an AFM^^. 

This manuscript is organized as follows. In Sect. [n| 
we present the model we have used to describe a pair of 
Coulomb-coupled MD2DESs, and introduce the linear- 
response functions which describe collective electron dy- 
namics. In Sect. |III| we present and discuss our main 



analytical and numerical results for the dispersion of op- 
tical and acoustic plasmons in these systems. Finally, in 



Sect. IV we present a summary of our main conclusions. 



II. MODEL HAMILTONIAN AND RANDOM 
PHASE APPROXIMATION 
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with {{A^B))^ the usual Kubo product. Within the 
RPA these functions satisfy the following matrix equa- 
tion. 



(4) 



where Xo{q^ o;) is a 2 x 2 diagonal matrix whose elements 
xf'\q^(^) are the well-knowiP^^^ noninteracting (Lind- 
hard) response functions of each layer at arbitrary dop- 
ing Ui. The off-diagonal (diagonal) elements of the ma- 
trix V = {Vw}i,e=i,2 represent inter-layer (intra-layer) 
Coulomb interactions. 

The bare intra- and inter-layer Coulomb interactions 
are influenced by the layered dielectric environment (see 
Fig. [l]) . A simple electrostatic calculation^^ implies that 
the Coulomb interaction in the i = 1 (top) layer is given 

by 
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[(e2 + e3)e«^ + (ea - £3)6" 



(5) 



where 



We consider two unhybridized MD2DESs separated by 
a finite distance d and embedded in the dielectric environ- 
ment depicted in Fig.[T] The two systems are assumed to 
be coupled solely by Coulomb interactions. The Hamil- 
tonian describing this system reads^^ {h = 1) 

k,£,a,f3 

^ JqYI ^^^'iQ)pQ,^P-Q,^' ' (1) 

Here v is the bare Dirac velocity, taken to be the same in 
the ^ = 1,2 tunnel- decoupled layers, S is the area of each 
layer, Vu'{q) is the matrix of bare Coulomb potentials, 
and 

Pq/ = Y^i-qAa^kAa (2) 

is the density-operator for the ^-th layer. The Greek 
letters are honeycomb-sublattice-pseudospin labels and 



D{q) = [(ei+e2)(e2 + e3)e^' + (ei-e2)(e2-e3)e-^^] . (6) 

The Coulomb interaction in the bottom layer, ^22(^)5 can 
be simply obtained from Vii{q) by interchanging 63 ^ ei. 
Finally, the inter-layer Coulomb interaction is given by 

Vu{q) = V2i{q) = ^^e2. (7) 

Notice that in the "uniform" ei = 62 = 63 = e limit 
we recover the familiar expressions Vii{q) = V22{q) 
2^6^ /{eq) and VM = ^2i(g) ^ Fii(g) exp(-gd). Pre- 
vious work on TI thin film and DLG collective modes has 
assumed this limit, which rarely applies experimentally. 



HI. COLLECTIVE MODES 

The collective modes of the system described by the 
model Hamiltonian ([T]) can be determined by locating the 
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poles ofY(g,cj) in Eq. Q. A straight forwa rd inversion 
of Eq. Q yields the following conditior]^^^ 



e{q,u;) = [l-Vi^{q)xf\q,oj)][l-V22{q)x'i'\q,oj)] 



(0), 



(8) 



The collective modes occur above the intra-band particle- 
hole continuum where x^^^ is real, positive, and a decreas- 
ing function of frequency. Eq. (|8|) admits two solutions. 



a higher frequency solution^ 



at a;op(^) which corre- 



sponds to in-phase oscillations of the densities in the two 
layers, and a lower frequency solution at Co'ac(^) which 
corresponds to out-of-phase oscillations. 

The plasmon collective modes of MD2DESs are of spe- 
cial interest because of the ease with which they may be 
altered by changing the carrier densities in either layer 
using gates. We note in particular that the carrier densi- 
ties in different layers can easily differ radically. For this 
reason we present our results in terms of the total 2D 
carrier density n = ni +n2, and the density polarization 
C = (^2 — ni)/n G [—1, 1]: C = 1 when the carrier den- 
sity is non-zero only in the bottom layer (ni = 0), while 
C = when the two layers have identical carrier densities 
(ni 712). 



A. Analytical results 

In this Section we report on exact analytical expres- 
sions for the RPA optical and acoustic plasmon dis- 
persions co'op,ac(^) that are valid in the long- wavelength 
g ^ limit where co'op(<7 ^ 0) cx and co'ac(<7 ^ 0) cx 

We start by deriving an exact expression for the RPA 
long-wavelength acoustic-plasmon group velocity. 



Cs = lim 

q^O q 



(9) 



Following Santoro and Giuliani^^, we first introduce the 
power expansion 



(10) 



for the acoustic-plasmon dispersion relation, and then 
define a function 



F{q) = e{q, c^q + Cag^ + c^q^ + ...)• 



(11) 



In the limit q ^ {) the function F{q) has the following 
Laurent-Taylor expansion 



F{q) = f-i q-' + fo + fiq + f2q'' 



(12) 



where the coefficients fi can be extracted from the ana- 
lytical expressioiP^M^ for the MD2DES Lindhard func- 
tion xf\Q^^)' For Eq. ^ to be valid we have to require 
that the coefficients fi vanish identically. The coefficient 
/_i depends only on Cg and by equating its expression 
to zero we arrive after some tedious but straightforward 



algebra at the following equation for x = Cs/v, the ra- 
tio between the plasmon group velocity Cg and the Dirac 
velocity v: 



2gsg.a,J{C^ -!)[! + 2x{Vx^-l - x)] 
- V2e2[l + x{Vx^ - 1 - x)]f{C) = , (13) 

where (g^) are real-spin (valley) degeneracy factors. 
In the case of DLG, g^ = g^ — 2^ while in the case of 
thin- film TIs = = 1- Iii Eq. (13) d = dk^ is a 



di mensionless inter-layer distance calculated with 
^/4:7^n/{gsgv) and n = ni + n2, and 



/(c) = (i+c)\/w + (i-c)\/rTc. (14) 

Eq. (p^Sl) can be conveniently solved for x by making the 



change of variables x T = \Jx^ — \ — x. After some 
straigthforward algebra we find that 



Cg 1 + A(aeeC^/e2,C) 

V ~ [l + 2A(aeeJ/e2,C)]^/' 



with 



A(aeeC^/e2,C) 



^g^v ^2(1 - O (^eed 



£2 



(15) 



(16) 



Eqs. (15)-(16) are the principle results of this Article. 
We see from this analytic expression that Cg is indepen- 
dent of ei and 63 and depends only on the barrier mate- 
rial dielectric constant, which in the case of TI thin films 
is simply the TI bulk dielectric constant. This behav- 
ior is a consequence of the out-of-phase character of this 
collective mode in which the double-layer total charge 
is locally constant but shifts dynamically between lay- 
ers. Because TIs tend to have narrow gaps they tend to 
have large dielectric constants (e2 ~ 100 in the case^^ of 
Bi2Te3). Thin- film collective modes will therefore tend 
to have Cq/v values that are quite close to 1 unless d is 
very large. (For large d the long- wavelength limit for- 
mula, which applies when both qd and g/Zcp are small, 
will have a limited range of applicability.) 

It follows from Eq. (15) that the ratio Cs/v is larger 
than unity for any value of the parameters o^ee, d^ 
and 62 so that the acoustic plasmon always lies out- 
side of the MD2DES particle-hole continuum. This im- 
plies than the acoustic plasmon is strictly speaking never 
Landau damped at small q. (A similar conclusion was 
reached previousl}^^ for the case of conventional 2D elec- 
tron gases, but was limited to the case of identical density 
and hence identical Fermi velocity.) 



For moderate values of d, however, Eq. (15) predicts 



a TI thin film sound velocity so close to the top of the 
particle-hole continuum that it will likely be unobserv- 
able because of damping effects not captured by the RPA, 
and because of disorder, which is always present to some 
degree. For the case of DLG, on the other hand, we ex- 
pect that acoustic plasmon collective modes will be well 
defined. This is particularly true in the case of DLG 
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FIG. 2: (Color online) Long- wavelength acoustic plasmon dis- 
persion of Coulomb- coupled massless-Dirac two-dimensional 
electron systems. The circles and squares are acoustic plas- 
mon frequencies oJadq) (in units of Sf = vkp) as functions 
of q/kF calculated numerically from the solution of Eq. ([8|. 
Here /cf is the Fermi wav e vector evaluated at the total density 
n = 721+712, i.e. /cf = Y^47r(ni + n2) / (gsdv) ■ The parameters 
we have used to calculate the curve labeled by — 0" ^^-re: 
^3 = = 2, m = 712 = 5 X 10^^ cm~^, ttee 2.2, d 3.35 A, 
ei = 62 = 1, and 63 = 3.9. These parameter values correspond 
to the case of two graphene layers on Si02 that are decoupled 
by rotation. The data labeled by = 0.8" have been calcu- 
lated by setting 721 = 1 x 10^^ cm~^ and 722 = 9 x 10^^ cm~^ 
with the same values for the other parameters. The solid lines 
plot uj — Csq for the C — ^ ^ind = 0-8 cases with the plas- 
mon group velocity Cs calculated from the analytical result, 
Eq. (15). These numerical results confirm the validity of our 



analytic result for Cs and the importance of accounting for the 
delicate dependence of the long- wavelength Lindhard function 
onv^uj/ (vq). The dashed line plots the upper-bound of the 
intra-band electron-hole continuum, cj = vq. 



with a small number of layers of BN as barrier mate- 
rial. When the BN barrier layer is very thin, the use 
of macroscopic dielectric parameters to characterize its 
screening properties is approximate; in that case mea- 
surement of the acoustic plasmon group velocity com- 
bined with Eqs. (p^Sll-iplGl) would allow the effective value 



of 62 to be determined experimentally. 

We note that an analytic result for Cg was reported 
previously in Ref. 20 [see their Eq. (5b)] for the special 
case of DLG embedded in a uniform dielectric, i.e. for 
ei = 62 = 63. In our notation, their result reads 



HDS 



272^1^ 



.-,1/2 



vT^ + vTTC ^2 



(17) 



This equation is evidently different from Eq. (15) above. 



We believe that Eq. (15) is the correct RPA result for 



the acoustic-plasmon group velocity and that Eq. (17) is 



of the Lindhard function xf\Qi^) as a function of wave 
vector q and frequency u in the region in which both 
these quantities are small. (See Sect. 4.4.3 of Ref. [T6l ) 

In particular, the limit of xf^\Q^^) ioi q ^ and u ^ 
depends on the ratio u = uo/{vq)^ i.e. on the direction 
along which the origin of the (g, cj) plane is approached: 
different limits are obtained for different values oi v. In 
an acoustic plasmon, the ratio v approaches a constant 
as g ^ and thus the limit of X^i\(l-> ^) which matters is 
the one in which g ^ while the ratio u/q is kept con- 
stant. This is the limit we have takerP^in the derivation 
of Eq. (pis]) - see Eq. ( pTj ). Eq. (17) is obtained by incor- 



rectly letting q ^ while uj is kept constant [see Eq. (4) 
in Ref. 20 : in this limit u diverges instead of going to a 
constant value. 

A careful comparison between our analytical prediction 
in Eq. (15) and the result obtained by the brute-force 
numerical solution of Eq. (|8| is shown in Fig. |2] We 
clearly see that Eq. ( 15 ) compares very well with the full 
numerical result. 

The analytical analysis of the long-wavelength opti- 
cal plasmon mode is simpler since this mode satisfies 
0Jop{q) ^ -s/Q for g ^ and therefore occurs at = 
uo/ivq) 00. We obtain an analytic result using the 
well-known high-frequency {00 ^ vq and uj <^ 2£f/) dy- 
namical limit of Xi^\Qi^)'- 



V (0)/ \ q 
hm Xi ^q^uj) = gs9v-. 2 



(18) 



with sy^i = vk^^i = v^ATTni/{g^g^) 
Eq. (Isl) we immediately find 



Using Eq. (18) in 



"op 



(^^0) 



V / 

2e 




with e = (ei +63)72. Note that Eq. ([19| does not depend 
on the inter-layer distance or on the dielectric constant 
62 , but only on the average e between top and bottom di- 
electric constants. Notice also that, in the limit ni ^ 
(C = 1), Eq. (19) reduces to the well-k nown p lasmon fre- 



incorrect. The difference is due to the singular behavior 



quency in a single-layer graphene sheet^^SlMj with elec- 
tron density 77,2. This expression applies for qd 1, in 
which case the entire double-layer MD2DES acts in the 
optical plasmon mode like a single conducting layer at 
the interface between dielectric media characterized by 
constants ei and 63. 



B. Numerical results 

In this Section we briefly report some representative 
numerical results for the optical and acoustic plasmon 
dispersion relations obtained by solving Eq. ([8|, dis- 
cussing first DLG and then TI thin films. 

In Fig. [3] we illustrate the typical properties of DLG 
plasmon modes for the case with the smallest MD2DES 
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FIG. 3: (Color online) Panel a) Optical and acoustic plas- 
mon dispersions (in units of the Fermi energy sy = vk^) 
in a twisted double-layer graphene system on a Si02 sub- 
strate as functions of wave vector q [in units of /cf = 
yjAn(ni + n^) / {g^gv)] • The values of the parameters that we 
have used to produce the data in this figure are: g^ — g^ — 2^ 
ni = 722 = 5 X 10"^^ cm~^ (corresponding to n = 10"^^ cm~^ 
and C, — 0), a^e — 2.2, d — 3.35 A, ei = 62 = 1, and 63 — 3.9. 
The intersections between the plasmon dispersions and the 
short-dashed line give the critical wave vector qc at which 
Landau damping starts. Panel b) Same as in panel a) but 
for C, = 0.5 [in producing the data shown in panel b) we have 
fixed the total density at the value used to produce the data 
in panel a), i.e. n = 10^^ cm~^]. 



separation, two adjacent layers on a Si02 substrate (ei = 
62 = 1 and 63 = 3.9) that are weakly hybridized e.Q^_ be- 
cause of a twist between their orientations^^. Fig. [3p) is 
for a symmetric system with the same electron concen- 
tration on the two layers (^ = 0), while Fig. ^p) refers 
to a system with a 50% density imbalance. The char- 
acteristic behaviors cJop(^) cx y/q of the optical plasmon 
and co'ac(^) q oi the acoustic plasmon are clearly visi- 
ble. The collective modes are not Landau damped when 



they appear in the gap between intra-band and inter- 
band particle-hole continua. When the two layers have 
different densities, their particle-hole continua are differ- 
ent and the gap is smaller for the lower density layer. For 
adjacent but twisted DLG systems d is small even when 
the carrier density is large {d ^ 0.2 in Fig. [3]). It follows 
that qd is small and the two MD2DESs are strongly cou- 
pled over the entire relevant frequency regime. In this 
small d example the acoustic plasmon frequency is close 
to the particle-hole continuum because the capacitive en- 
ergy associated with charge sloshing between the layers 
is proportional to the small layer separation. 

In Fig. [4] we illustrate the strength of plasmon decay by 
emission of single electron- hole pairs (Landau damping). 
Notice that Landau damping occurs when the curves 
^op,ac(<?) in Fig.|3]hit the inter-band electron-hole contin- 
uum of the layer with lower density (layer "1" in our con- 
vention). The larger the sooner this happens. In par- 
ticular, in the limit in which layer "1" is neutral (^ = 1), 
Landau damping is present from vanishingly small wave 
vectors: damping of the optical plasmon excitation asso- 
ciated with electrons in the high-density layer starts at 
arbitrarily small wave vectors since decay can easily occur 
via the emission of inter-band electron-hole pairs in the 
neutral layer. The many-body properties of two or more 
decoupled graphene layers can thus be strongly affected 
by inter-layer Coulomb interactions, even by apparently 
innocuous geometric features such as the presence of a 
nearly-neutral layer. 

In Fig. [5] we compare optical and acoustic plasmon 
dispersions for DLG and TI thin-film systems. For the 
TI thin-film case we have chosen the following param- 
eters: i) ei = 1, 62 = 100 (this roughly corresponds 
to the dielectric constant of Bi2Te3), and 63 = 4.0; ii) 
a total electron density on the top and bottom surface 
states of n = 10^^ cm~^; and iii) a thickness of the 
TI slab of = 6 nm, corresponding to a six quintile 
layer MBE-grown Bi2Te3 film (J ^ 6.7). The DLG ex- 
ample has the same total density and layer separation 
(J ~ 3.4; the difference in d in the two cases stems from 
the ^s/^v spin/ valley degeneracy factors) and dielectric 
constants ei = 1 and 62 = 63 = 4.0, corresponding to 
two graphene layers separated by approximately 15 BN 
layers and lying on a BN substrate. In both cases we see 
that a crossover occurs at intermediate values of q be- 
tween strong (small q) and weak (large q) coupling of the 
two collective modes. In the TI case the higher frequency 
optical plasmon mode deviates much more strongly from 
simple -s/q behavior at this crossover because strong di- 
electric screening by the TI bulk suppresses the single- 
surface plasmon mode. [Note however that the effective 
dielectric constant for this limit is (e2 + ei^3)/2 rather 
than 62 as used in Ref. |2T1] The acoustic plasmon mode 
of the TI thin film case is, on the other hand, strongly 
suppressed in the strong-coupling limit, as discusses ear- 
lier, and has a velocity much closer to the bare Dirac 
velocity than in the corresponding DLG case. 
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FIG. 4: (Color online) Landau damping of collective modes 
in double-layer graphene. Panel a) The absolute value of the 
imaginary part of the Lindhard function of the top layer, 
l^m Xi^\q^^)\j evaluated at the frequency uo = ^op(^) 
[lu = U0ac{q)] of the optical [acoustic] plasmon. The data in 
this plot refer exactly to the parameters used in Fig. [3^). 
Note that, within RPA, Xi^'* (^5 '^op,ac(^)) is identically 
zero for wave vectors ^ up to a critical value gop,ac at which 
<^op,ac(^) hits the inter-band electron-hole continuum asso- 
ciated with the low-density layer. Since data in this panel 
correspond to ^ = O5 top-layer and bottom-layer Lindhard 
functions are identical. Panel b) Same as in panel a) but for 
( = 0.8 (m = 1 X 10^^ cm"^ and 712 = 9 x 10^^ cm"^). Note 
that the ^op,ac decreases with increasing becoming zero in 
the limit ( ^ 1. Since in this panel C / we have plotted 
both top-layer (low-density) and bottom-layer (high-density) 
Lindhard functions. 



FIG. 5: (Color online) Optical and acoustic plasmon disper- 
sions (in units of the Fermi energy sf = vkp) in a double- 
layer graphene system [panel a)] and a topological insulator 
thin- film [panel b)] as functions of wave vector q [in units of 
kp = i/47r(ni + n2)/ (gsgv)]- The intersections between the 
plasmon dispersions and the short-dashed line give the criti- 
cal wave vector Qc at which Landau damping starts. Panel a) 
The values of the parameters that we have used to produce the 
data in this figure are: = 5'v = 2, ni = 712 = 5 x 10"^^ cm~^ 
(corresponding to n = 10"^^ cm~^ and C = 0), <^ee = 2.2, 
d = 6 nm, ei = 1, and £2 = £3 = 4.0. Panel b) The val- 
ues of the parameters that we have used to produce the data 
in this figure are: = 5'v = 1, '/T'I = 712 = 5 x lO"*^^ cm~^, 
aee — 4.4, d = 6 nm, ei = 1, £2 = 100, and €3 — 4.0. Note 
that due to the large value of the bulk TI dielectric constant, 
the acoustic plasmon is almost locked to the top of the intra- 
band electron-hole continuum. 



IV. DISCUSSION AND CONCLUSIONS 

We have presented an analysis of the electronic collec- 
tive modes of systems composed of two unhybridized but 
Coulomb-coupled massless-Dirac two-dimensional elec- 
tron systems (MD2DESs) separated by a vertical dis- 



tance d. The primary example we have in mind is topo- 
logical insulator (TI) thin films, which are always de- 
scribed at low energies by this type of model because 
topologically protected MD2DESs always appear on both 
top and bottom surfaces. Also of interest are closely re- 
lated systems, which we refer to as double-layer graphene 
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(DLG) systems, containing two graphene layers that are 
weakly hybridized either because they are rotated rela- 
tive to each other or because they are separated by a 
dielectric barrier layer. Importantly, we allow for a gen- 
eral dielectric environment in which the material above 
the top MD2DES layer (ei), between the two layers (62), 
and below the bottom MD2DES layer (63) are all allowed 
to have different dielectric constants. In the case of TI 
thin film 62 is the bulk dielectric constant of the TI which 
is expected to have large values. The carrier collective 
modes of MD2DESs are expected to be most robust in 
the gap between intraband and interband particle-hole 
excitations. 

The double-layer systems of interest quite generally 
have two collective modes which in the limit of small 
qd involve density-fluctuations in the two-layers that are 
strongly coupled, and in the limit of large qd weakly cou- 
pled single-layer plasmons. One key parameter which 
controls collective mode properties is the dimensionless 
product k^d = d. Small values of k^d imply that the 
layer separation is smaller than the typical distance be- 
tween electrons within a layer and that collective modes 
at all values of up to ^ kp are strongly coupled com- 
binations of the two individual layer density-fluctuation 
contributions. For large k^d a. crossover occurs for q G 
(0, /cp) between strongly and weakly coupled double-layer 
collective modes. Both small and large values of d are 
achievable in samples where disorder plays an inessential 
role in both DLG and TI thin film cases. 

Our study focuses on the long-wavelength limit in 
which both qd and q/k^ = qd/d are small. We have 
derived analytic expressions for both frequencies of both 
the low-energy linearly dispersing acoustic plasmon mode 
^ac(^) and for the high-energy optical plasmon mode 
co'op(^) which has ^ dispersion at long- wavelengths. In 
this limit we find that uJac{Q) — vq ex l/e2 whereas 
'^op(^) Y^2/(ei + 63); i.e. the separation of the acous- 
tic plasmon mode from the upper edge of the intra-band 
particle-hole continuum is very strongly suppressed by 
a large bulk TI dielectric constant, whereas the coupled 
double-layer plasmon mode is unaffected. This double- 
layer optical plasmon behavior contrasts with that of 
a large qd single-surface plasmon mode which has a 
frequency proportional to y^2/{e2 + ^1^3). The long- 



wavelength limit of cjac {q) is sensitive not only to the en- 
ergy associated with inter-layer charge sloshing but also 
to its microscopic kinetics as captured by the singular 
sensitivity of the MD2DES Lindhard function to uj/{vq). 
By carefully accounting for this dependence we are able 
to correct a previous analytic expression in a way that is 
quantitatively particularly important in the TI thin film 
(large 62) case. 

Double-layer collective mode coupling plays an im- 
portant role in MD2DES correlations when d is small. 
Even when J is large, strongly-coupled small qd = qd/k^ 
modes will often be experimentally accessible and may 
play an important role in graphene multi-layer or TI 
based plasmonics. The analytic results derived in this 
paper can be used to readily anticipate how these modes 
depend on system parameters. 

From the more theoretical point of view, it will be 
intriguing to study physical properties of plasmons in 
Coulomb-coupled MD2DESs beyond the random phase 
approximation by employing e.g. many-body diagram- 
matic perturbation theor}!^. 
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